REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis 
Highway,  Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a 
collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1 .  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

17-05-2011 

3.  DATES  COVERED  (From  -  To) 

4.  TITLE  AND  SUBTITLE 

A  Generic  4th  Order  2D  Unstructured  Euler  Solver  for  the  CESE  Method 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

David  L.  Bilyeu,  Yung-Yu  Chen  and  S.-T.  John  Yu 

5d.  PROJECT  NUMBER 

5f.  WORK  UNIT  NUMBER 

23041057 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory  (AFMC) 

AFRL/RZSS 

1  Ara  Road 

Edwards  AFB  CA  93524-7013 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

AFRL-RZ-ED-AB-201 1-167 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory  (AFMC) 

AFRL/RZS 

5  Pollux  Drive 

Edwards  AFB  CA  93524-7048 

10.  SPONSOR/MONITOR’S 
ACRONYM(S) 

11.  SPONSOR/MONITOR’S 
NUMBER(S) 

AFRL-RZ-ED-AB-201 1-167 

12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited  (PA  #11 129). 

13.  SUPPLEMENTARY  NOTES 

For  the  AIAA  Aerospace  Sciences  Meeting,  Nashville,  TN,  9-12  Jan  2012. 

14.  ABSTRACT 

Previously,  Chang  reported  a  new  class  of  high-order  Conservation  Element  Solution  Element,  CESE,  methods  for  solving 
nonlinear  hyperbolic  partial  differential  equations.  The  scheme  was  then  extended  by  Bilyeu,  et  al.  to  solve  a  ID  vector 
equation  with  an  arbitrary  order  of  convergence.  This  new  high-order  CESE  method  shares  many  favorable  attributes  of 
the  original  second-order  CESE  method,  including:  (i]  compact  mesh  stencil  involving  only  the  immediate  mesh  nodes 
surrounding  the  node  where  the  solution  is  sought,  (ii]  the  CFL  stability  constraint  remains  the  same,  i.e.,  _  1,  as ,  and  [iii] 
superb  shock  capturing  capability  without  using  an  approximate  Riemann  solver.  In  the  present  extended  abstract,  we 
extend  the  ID  formulation  to  2D.  A  general  formulation  is  presented  for  solving  the  coupled  equations  with  4thorder 
accuracy.  To  demonstrate  the  formulation,  two  linear  cases  are  reported.  The  linear  test  cases  shown  are  the  convection 
equation  and  the  linear  acoustic  equation.  In  the  final  paper  we  plan  to  present  non-linear  test  cases  and  solve  the 
equations  on  a  quad  mesh. 


15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT 

18.  NUMBER 
OF  PAGES 

19a.  NAME  OF  RESPONSIBLE 
PERSON 

Dr.  Jean-Luc  J.  Cambier 

a.  REPORT 

Unclassified 

b.  ABSTRACT 

Unclassified 

c.  THIS  PAGE 

Unclassified 

SAR 

13 

19b.  TELEPHONE  NUMBER 

(include  area  code) 

N/A 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  239.18 


Unclassified 


Unclassified 


Unclassified 


A  Generic  4thOrder  2D  Unstructured  Euler  Solver  for 

the  CESE  Method 

David  L.  Bilyeu,*  Yung- Yu  Chen,  t  and  S.-T.  John  Yu  * 

The  Ohio  State  University,  Columbus,  OH  43210,  USA 


Previously,  Chang  reported  a  new  class  of  high-order  Conservation  Element  Solution 
Element,  CESE,  methods  for  solving  nonlinear  hyperbolic  partial  differential  equations. 
The  scheme  was  then  extended  by  Bilyeu,  et  al.  to  solve  a  ID  vector  equation  with  an 
arbitrary  order  of  convergence.  This  new  high-order  CESE  method  shares  many  favorable 
attributes  of  the  original  second-order  CESE  method,  including:  (i)  compact  mesh  stencil 
involving  only  the  immediate  mesh  nodes  surrounding  the  node  where  the  solution  is 
sought,  (ii)  the  CFL  stability  constraint  remains  the  same,  i.e.,  <  1,  as  ,  and  (iii)  superb 
shock  capturing  capability  without  using  an  approximate  Riemann  solver.  In  the  present 
extended  abstract,  we  extend  the  ID  formulation  to  2D.  A  general  formulation  is  presented 
for  solving  the  coupled  equations  with  4thorder  accuracy.  To  demonstrate  the  formulation, 
two  linear  cases  are  reported.  The  linear  test  cases  shown  are  the  convection  equation  and 
the  linear  acoustic  equation. 

In  the  final  paper  we  plan  to  present  non-linear  test  cases  and  solve  the  equations  on  a 
quad  mesh. 


I.  Introduction 

In  this  work,  we  extend  the  ID  solver  for  one  nonlinear  hyperbolic  equation  to  two  dimensions.  The 
new  formulation  is  currently  restricted  to  4thorder  but  can  be  extended  to  higher  orders  by  extending  the 
Taylor  series.  To  demonstrate  the  capabilities  of  the  new  scheme,  we  apply  the  method  to  solve  two  sets  of 
equations:  (i)  the  linearized  acoustic  equations,  and  (ii)  a  convection  equation. 

The  original  second-order  CESE  method  ’  1  solves  the  hyperbolic  PDEs  by  discretizing  the  space-time 
domain  by  using  the  conservation  elements  (CEs)  and  solution  elements  (SEs).  The  profiles  of  unknowns  are 
prescribed  by  assumed  discretization  inside  SEs.  Aided  by  the  approximation  for  the  unknowns  in  the  SEs, 
space-time  flux  conservation  can  be  enforced  over  each  CE.  The  calculation  of  space-time  flux  conservation 
results  in  the  formulation  for  updating  the  unknowns  in  the  time  marching  process.  The  special  features 
of  the  CESE  method  include:  (i)  The  conserved  variables  are  represented  by  a  Taylor  series  expansion  in 
both  space  and  time.  The  order  of  the  Taylor  series  is  also  the  order  of  the  solver,  (ii)  The  a  scheme,  the 
core  scheme  of  the  CESE  method,  is  non-dissipative.  (iii)  The  CESE  method  has  the  most  compact  mesh 
stencil  possible  involving  only  the  immediate  neighboring  mesh  points  that  surround  the  mesh  node  where 
the  solution  is  sought,  (iv)  The  method  uses  explicit  integration  in  time  marching.  The  stability  criterion  is 
CFL  <  1.  (v)  No  approximate  Riemann  solver  is  used  and  the  scheme  is  simple  and  efficient. 

The  space-time  stencil  used  in  this  derivation  is  the  same  as  that  reported  by  Wang  and  Chang J  and  is 
repeated  here  for  completeness.  Figure  1  is  a  top  down  view  of  the  computational  domain.  In  this  figure  the 
solid  lines  represent  the  mesh  and  the  dashed  lines  are  the  boundaries  of  the  conservation  element  centered 
at  G(j).  The  squares  are  vertices’s  of  the  mesh,  the  circles  are  the  solution  points  and  the  x  are  the  centroid 
of  each  triangle.  Figures  2a  and  2b  are  the  conservation  element,  CE,  and  solution  element,  SE,  associated 
with  the  solution  point  G(j).  These  CE  is  the  3D  space-time  element  that  is  used  to  conserve  flux  and  the 
SE  is  the  domain  that  that  the  Taylor  series  is  assumed  to  be  valid  in. 

In  the  following  derivation  we  make  the  following  assumptions: 
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tPh.D.  Candidate,  Dept,  of  Mechanical  Engineering,  chen.1352@osu.edu. 
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1.  The  flux  is  a  known  function  of  the  flow  variables 

2.  Inside  of  a  solution  element  the  flow  variables  and  fluxes  can  be  expressed  as  a  Taylor  series 

3.  The  desired  order  of  convergences  is  even 

4.  Each  even  derivative  has  an  odd  derivative  in  each  spatial  direction. 

5. 


d2u  d2u 
dxdy  dydx 


6.  The  source  term  is  handled  via  operator  splitter 

7.  For  non  steady  state  problems  the  order  of  convergence  of  the  initial  conditions  must  match  or  exceed 
the  desired  order  of  accuracy  of  the  solver. 

Shown  below  is  a  4thorder  Taylor  expansion  in  two  spatial  derivatives 


u*  =  u  +  ux Ax  +  uvAy  +  (uxx Ax2  +  uxyAxAy  +  uyxAyAx  +  uvyAy2}  /2.0+ 
( ' V'xxxAx  T  UyyyAy  -f-  ( uxxy  T  uxyx  T  'Uyxx') Ax  Ay  T  ( uxyy  T  uyxy  T  'Uyyx') Ay  AxJ  /6.O. 

If  we  assume  that  the  Taylor  series  expansion  of  uxy  and  uyx  are  equal  then 

duXy  diiyx  ,  duxy  duyx 

dx  dx  a  dy  dy 

or 

UXyx  —  'Uyxx  and  UXyy  —  Uyxy. 


(1) 


(2) 


Under  theses  assumption  the  Taylor  series  becomes 

u*  =  u  +  ux Ax  +  uyAy  +  (uxx Ax2  +  2uxyAxAy  +  uvvAy 2)  /2.0+ 

(j^xxxAx  4“  UyyyAy  T  ( uxxy  ~f“  2 uxyx) Ax  Ay  -t-  (2 uxyy  T  uyyx) Ay  Ax )  /6.O. 

This  extended  abstract  is  organized  as  follows.  Section  II  derives  the  equations  necessary  to  progress  the 
conserved  variables  to  the  next  time  step.  This  section  is  broken  up  into  4  sub  sections  one  for  calculating 
the  fluxes,  temporal  derivatives,  even  and  odd  derivatives  of  the  conserved  variables.  Section  III  contains 
the  preliminary  results  from  the  convection  equation  and  linear  acoustic  equation. 


II.  Derivation 


The  governing  equation  that  we  wish  to  solve 


is 


dU 

dt 


-  v  .  -pXj  =  S. 


where 


(4) 


U  =  (Ul,  «2, . . . ,  um)T,  Fx  =  (/f ,  /£, . . . ,  /£f,  Fy  =  C ff,  fl  ■  ■  • ,  /^)T,  S  =  (Sl,  s2, . . . ,  sm)T , 

m  is  the  number  of  equations.  Each  Um  is  expressed  as  a  Taylor  series, 

N  N-aN-a-b 


-EE  E 

a= 0  6=0  c= 0 


Qa+b+c 


l 


dxadybdtc  a!6!c! 


(x-Xj)a(y-yj)b{t-tn)c, 


(5) 

(6) 


where  N  is  equal  to  the  desired  order  of  convergence  minus  1.  For  4thorder  N  is  equal  to  3.  For  future 
reference  the  power  in  the  numerator  of  the  d  will  be  dropped,  e.g. 


da+b+cum  =  dum 
dxadybdtc  dxadybdtc 


(7) 
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The  derivative  of  any  Taylor  series  can  be  written  as 


where  A  = 
is  also  true 


du* 


dxIdyJdtK 


A  B—a C—a—b 

=  EE  E 


dur 


1 


dxI+a8yJ+bdtK+c  a\b\c\ 


(x-xj)a(y-yj)b{t-tn)c, 


(8) 


a= 0  6—0  c— 0 

N  —  I  —  J  —  K ,  B  =  N  —  I  —  J  —  K  and  C  =  N  —  /  —  J  ~  K.  A  similar  expansion  for  the  fluxes 


df*m 


dxIdyJdtK 


A  B—a C—a—b 

-EE  E 

a— 0  b— 0  c— 0 


dfrr 


1 


dxI+adyJ+bdtK+c  a\b\c\ 


(x  —  Xj)a(y  —  yj)b(t  —  tn)c 


(9) 


For  2D  there  are  a  total  of  3  o  (-^  —  n)(n  +  1)  unknowns  per  governing  equation.  For  a  4thorder  2D  Euler 
equation  there  would  be  a  total  of  180  unknowns.  In  what  follows  we  will  show  that  the  only  independent 
variables  are  the  spatial  derivatives  of  the  conserved  variables 


dllm 

dxIdyJ 


1  =  0,1,.. .,1V 
J  =  0, 1, . . . ,  IV  —  J 


to  =  1, ... ,  Neq. 


For  a  4thsolver  the  independent  variables  are 


(10) 


Table  1:  Table  of  the  independent  conserved  variables 


Even  Odd 


u 

ux 

Uy 

^ XX 

V'XXX 

'U'xxy 

'U'xy 

'U'xyx 

^ xyy 

uyy 

'Uyyx 

uyyy 

II. A.  Fluxes 


Since  the  fluxes  are  known  functions  of  the  conserved  variables  the  derivatives  of  the  fluxes  can  be  found 
from  the  chain  rule. 

Neq 

f  I  T  ~  /till 

Y1=x,y,t  (11) 


\  ' 


dfxm  dm 


dYi 


dm  dYi  ’ 


a2/r 

SYi<9y2 


Neq 

=  E 


dfm  d2m 
dm  dY\  dY 2 


duk 


“  duiduk  dY1  dY2 

L.K, 


(Yi,Y2)  = 


(x,x) 

(t,t) 

(x,t) 


(y,y) 

(x,y) 

(y,t) 


(12) 


8Y1dY2dY3 


Neq 

=E 

i 

Neq 


dft 


d3m 


dm  8YidY2dY3 
d2fm  (  d2m  duk 


“  dmduk 


Neq 

E- 

i,k,p 


d2m  duk  d2m  du/, 


d3K 


dY1dY28Y3 
dm  duk  dup 


dmdukdup  dY\  dY2  dY3 


dY]  8Y3  dY2  8Y2dY3dY1 


(x,x,x)  {y,y,y)  (t,t,t) 

(Y1,Y2,Y3)  =  (x,x,y)  (y,y,t)  (x,x,t) 

{x,y,y)  (y,t,t)  {x,t,t) 


(x,y,t) 


(13) 


for  xl  =  x,y.  These  equations  will  continue  until  all  of  the  fluxes  are  found.  After  these  equations  are 
applied  there  will  be  n  (IV  —  n)(n  +  1)  unknowns  left  per  governing  equation.  This  method  will  work 
for  all  higher  derivatives  but  becomes  more  and  more  complicated  as  the  number  of  derivatives  increases. 
Another  approach  is  to  directly  calculate  each  derivative  with  out  using  a  Jacobian. 
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II. B.  Temporal  derivatives 

In  this  section  we  will  use  the  governing  equation  and  its  derivatives  to  find  the  temporal  derivatives  of  the 
conserved  variables. 


dlL-m 

gS 

dfvm 

dt 

S-rri 

dx 

dy 

9  ( 

/  durn  \ 

dSm 

d2fxm 

d2Pm 

dx  \ 

v  dt  J 

dx 

dxdx 

dydx' 

d  / 

'dum\ 

dSm 

d2Pm 

d2Pm 

dy  \ 

,  dt  J 

dy 

dxdy 

dydy  ’ 

In  general  any  temporal  derivative  of  a  conserved  variable  can  be  found  using 

dum  dsm  df £  dfy, 


dxIdyJdtK  dxIdyJ  dtK~ 1 
for  K  =  1,2,..., N 


dxI+1dyJ  dtK 


-l 


■7  =  0,1,.. 


,N-K-  1  =  0,1,. 


dx1 dyJ+1dtK~1 ’ 
N  -K  -  J. 


(14) 


(15) 


This  leaves  us  with  ^2^=0(n  +  1)  unknowns  per  governing  equation.  For  a  4thorder  Euler  solver  this  would 
be  40  unknowns,  10  per  equation. 


II. C.  Space-Time  Flux  Conservation 

To  find  the  next  set  of  conserved  variables  we  conserve  flux  across  both  space  and  time.  It  should  be  noted 
that  the  geometry  of  the  Conservation  Element  and  Solution  Elements  are  the  same  as  the  ones  presented 
by  Wang  and  Chang.  By  using  the  Gauss-divergence  theorem  we  can  convert  the  governing  equation  into 
integral  form 


hm (x,y,t)  ■  ds 


smdV. 


(16) 


This  equation  will  be  used  to  find  all  even  derivatives,  e.g.  u,  uxx,  uxy,  uyy,  uxxxx, . . . .  The  rest  of  the 
unknowns  will  be  found  using  a  central  difference  like  procedure.  Eq.  (16)  can  be  used  to  find  um  and  its 
derivative  will  be  used  to  find  the  other  even  derivatives 


/ 


dhm 

dxIyJ 


(■ x,y,t )  ■  ds 


dsm 

dxIyJ 


dV. 


(17) 


II. D.  Even 

The  even  derivatives  are  found  by  conserving  the  flux  through  the  CE.  To  make  this  calculation  more 
manageable  we  will  separate  this  calculation  into  three  steps,  (1)  the  flux  through  a  side  face,  (2)  the  flux 
through  the  bottom  face  and  (3)  the  flux  through  the  top  face.  The  flux  through  the  side  and  bottom  faces 
are  calculated  using  the  neighboring  solution  points  while  the  flux  through  the  top  face  is  calculated  from 
the  current  solution  point. 

To  begin  we  will  calculate  the  flux  through  the  side  faces.  Figure  3  shows  the  faces  associated  with  the 
solution  point  jr.  The  coordinate  axis  is  located  at  cell  center  and  the  vertex  ( xo,yo )  is  the  solution  point. 

First  we  will  shift  the  coordinate  axis  to  the  neighboring  solution  point.  For  future  reference  we  will 
denote  the  actual  location  with  capital  letters,  (X,  V),  and  relative  coordinates  with  lower  case  letters, 
(a :,y).  The  flux  can  now  be  expressed  as 

d/m  xayb(t —  t™-1/2)0 
^  ^  ^  dxadybdtc  a\b\c\ 

a  b  c  * 
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where  x  =  X  —  Xjr  and  y  =  Y  —  Yjr .  Next  we  apply  the  following  coordinate  transformation 


x  =  u(x  i  —  Xo)  +  Xq 
y  =  u(yi  -  y0)  +  y0 


t  = 


0  <  u  <  1, 


At 


-v  +  t' 


i-l/2 


0  <  V  <  1. 


(19) 


The  normal  vector  is  equal  to  N  =  [^r(yi  —  yo),  —  ^{x\  —  xq),  0].  Using  the  coordinate  transformation  Eq. 
(18)  can  be  rewritten  as 


/Ei*~EEE  d^m  (u(xi ~ x°^ + x°^a (u(yi ~ + y°^b 


a  b  c 

The  flux  through  a  side  face  is  equal  to 


dxadybdtc 


alblcl 


jj  {fm*Ni)r  dudv  = 
2  ABC  rl  rl 

EEEE 


df%_  \n  1/2  (u(  xi  -  x0)  +  x0)a  ( u(y j  -  y0)  +  y0)b 


Jo  \dxadybdtc  .  , 

i=l  a  b  c  •'O  JO  \  a  /  j. 

The  flux  through  the  faces  associated  with  solution  point  jr  is 

2  ABC 

( Fm)jr  =  E  E  E  E 


a\b\c\ 


Nrjdudv. 


(f)c+1 


i=  1  a  b  c 
r*  1 


dxadybdtc  I  (c  +  l)!a!6! 

(iV-)0i  J  {u(x\  -  x0)  +  xo)a  { u{y 1  -  yo)  +  yof  du+ 

(n*)  [  (u(x2  -  Xo)  +  x0)a  (u(y2  -  yo)  +  yo)b  du 

V  /  02  Jo 


(20) 


(21) 


(22) 


The  total  flux  through  the  side  faces  is  X)r=i  (Fm)jr- 

Next  we  will  calculate  the  flux  through  the  top  and  bottom  surface.  This  requires  us  to  integrate  over  an 
arbitrary  quadrilateral.  The  only  assumption  that  can  be  made  about  the  quadrilateral  is  that  it  is  concave. 
As  with  the  side  face  we  shift  the  centroid  to  the  current  solution  point.  An  arbitrary  bottom  face  is  shown 
in  Fig.  4.  In  this  figure  point  1  is  the  centroid  of  the  adjacent  CE  and  point  3  is  the  centroid  of  the  CE  and 
the  current  solution  point.  Points  2  and  4  are  the  shared  vertices’s  of  the  mesh.  To  integrate  the  flux  through 
the  bottom  surface  we  split  the  quadrilateral  into  two  triangles  and  then  provide  a  coordinate  transformation 
that  will  convert  the  triangles  into  right  triangles.  The  area  of  the  triangles  are 


Ai  =  \ 

x3  -xi  y3  -  yi 

Aii  =  \ 

X2  -  Xi  y2  -  yi 

1 

1 

x3  -X!  y3-  yi 

The  transformation  equation  used  for  region  /  is 

ui  =  ux( x  -  xi)  +  uv(y  -  yi)  vi  =  vx(x  -  xi)  +  vy(y  -  yi). 
To  express  x  and  y  as  functions  of  u  and  v  we  set  (M3,  V3)  =  (1,0)  {ua,  va)  =  (1,1) 

xi  =  u(x 3  -  Xi)  +  v{xa  —  *3)  +  %i 
yi  =  u(y3  -  yi)  +  v(yA  -  yo)  +  yi, 


(23) 


(24) 


(25) 


where  0  <  v  <  u  and  0  <  u  <  1.  The  normal  is  equal  to  IV  =  [0,  0,  2 A/].  This  process  is  repeated  for  region 
II  yielding 


Xii  =  u(x2  -  Xi)  +  v(x3  -  X2)  +  X\ 
yi!  =  u(y2  -  yi)  +  v(y3  -  y2)  +  yu 


(26) 
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where  0  <  v  <  u  and  0  <  u  <  1.  The  normal  is  equal  to  N  =  [0,  0,  2Ajj\.  Since  the  first  two  terms  in  the 
normal  are  zero  the  flux  through  the  bottom  surface  is 


{UmT, 


n— 1/2 


=  2 


1  pu  A  B 


0  JO 


EE 

a  b 


l 


u*jAj  +  u^jAjjdudv  = 
\  n— 1/2 

OUm  r 


alb!  \dxadyb 


A i  {u{x3  -  Xx)  +  v{x4  -  x3)  +  Xx)a  (u{y3  -  yi)  +  v(y4  -  y3 )  +  yx)  + 


Au  (u(x 2  -  x4)  +  v(x3  -  x2)  +  xx)a  ( u{y2  -  yx)  +  v(y3  -  y2)  +  yx) 


dudv. 


(27) 


The  total  flux  through  the  bottom  faces  is  XEi  (Um)™~1/"2 .  The  flux  through  the  top  face  is  found  in  a 
like  manner.  Since  the  bottom  and  top  face  have  the  exact  same  shape  we  use  the  same  diagram  for  both 
the  top  and  the  bottom.  The  difference  is  that  the  corresponding  solution  point  is  located  at  point  3.  Which 
means  that  (x3,y3)  =  (0,0).  Using  this  the  flux  through  the  top  face  is 


(Um)jr  —  2  /  /  u*jAj  +  u*jjAjjdudv  — 


1  nU  A  B 


0  */o 


EE' 

a  b 


dun 


alb !  \dxadyb 


Ai  ( u(x3  -  xx)  +  v(x4  -  £3)  +  xx)a  (u{y3  -  yx)  +  v(y4  -  y3 )  +  yx)  + 


An  ( u(x2  -  xx)  +  v(x3  -  x2)  +  Xx)a  ( u(y2  -  yx)  +  v{y3  -  y2)  +  yx) 


dudv. 


(28) 


The  total  flux  through  the  top  face  is  equal  to  X/=i  ■  These  equations  are  combined  to  progress  the 


even  derivatives  to  the  next  time  step 

3 


E  (U-V.  =  ///  s-dV  -  E  ((^r 1/2  +  w~>*)  • 


(29) 


Since  the  geometry  of  the  CE  and  SE  are  the  same  as  in  the  second  order  version  we  can  use  the  same 
simplifications.  This  results  in  the  first  derivatives  ux,uy  canceling  out.  This  allows  us  to  solve  for  um 
explicitly 


^  =  lb  (///  SmdV  -  E  ((<E)”  +  ( UmT-1' 2  +  {Fm)^  J  • 


Atot 

where  Atot  is  the  total  area  of  the  CE  and 

n:  - 


(30) 


0  ^0 


A  B 

EE 

a  b 


dur 


1 


alb!  \dxadyb 


Af  (u(x3  -  Xx)  +  v(x4  -  *3)  +  xx)a  ( u{y3  -  yx)  +  v(y4  -  y3)  +  yx)  + 


(0,0), 


An  (u(x2  -  xx) +  v(x3  -  x2)  +  xx)a  {u(y2  -  yx) +  v(y3  -  y2)  +  yx)  dudv(a,b)  ±  (1,0),  ■ 

(0,1) 


(31) 


None  of  the  integrals  in  Eqs  (22),  (27),  and  (28)  have  been  evaluated.  All  of  these  integrations  can  be 
found  using  multiple  methods.  The  easiest  method  is  to  use  a  program  capable  of  symbolic  integration. 
But  they  could  also  be  evaluated  by  hand  using  f  fdg  =  fg  —  f  gdf  multiple  times.  It  is  also  possible  to 
use  the  binomial  theorem  to  generate  a  formula  that  could  be  used  for  any  combination  of  a  and  b.  In  our 
investigation  we  evaluated  the  integrals  using  a  symbolic  math  package. 
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II. E.  Odd 


The  odd  derivatives  are  evaluated  using  a  central  difference  like  approach.  The  stencil  used  for  the  odd 
derivatives  is  shown  in  Fig.  5.  In  this  figure  the  black  dots  are  the  vertices’s  and  the  circles  are  the  solution 
points.  The  crosses  represent  the  shifted  locations  of  the  solutions  points.  A  shifted  point  is  used  to  reduce 
the  dissipation  as  the  CFL  number  approaches  zero.  The  location  of  the  crosses  is  calculated  by 

{Xj'r  >  Vj'r  )  =  [T(Xjr  ~  Xj)  +  Xj  >  T(Vjr  ~  Vj  )  +  Vj]  ,  (32) 

where  r  is  the  absolute  value  of  the  local  cfl  number.  We  take  a  Taylor  expansion  from  j  to  any  j' ,  r  =  1,  2,  3 
as 


(jV,n)  ££  (dxadyb)  .  a\b\  (Xjr  Xj )  (Yjr  Yj^  ' 

n= n h — n  x  »  /  7 


U„ 


a=0b=0  x~~  _i'  7  3 

We  will  shift  the  coordinate  axis  so  that  point  j  is  at  the  origin.  The  Taylor  expansion  now  becomes. 

A  B 


u. 


.ov.»>  =  ££(^) 

n—n h—n  \  &  /  7 


a= 0 6=0  v  ~  a  y  3 

In  order  to  find  the  two  unknowns  we  will  construct  two  equations. 


- 

Unix 

/  „  7 

(1) 

— 

X2  2/2. 

Umy 

(  9um 

n 

1  /6 

^  dxadyb 

1 

a\b\X  lV  1 

f  dum 

1  / a  fb 

y  dxadyb 

a\b\X  2  2/  2 

for  (a,  b)  ±  (0, 1),  (1,0). 


Since  can  not  be  evaluated  yet  we  approximate  it  with 
ABC 


a=0  6=0  c=0 


dur 


1 


/  At  V 


dxadybdtc  J  ■  alb\c\ 


( X'ir  -  Xir)  (■ Vjr  -  Vjr) 


(33) 


(34) 


(35) 


(36) 


Eqs.  (35)  and  (36)  are  used  to  find  one  possible  solution  to  umx  and  umy  using  information  from  two  of  the 
three  surrounding  solution  points.  This  will  be  repeated  using  all  possible  combinations  of  adjacent  solution 
points,  e.g.  (ji ,ji),  (J2J3),  (.7,3 Ji ) ■  These  solutions  are  then  weighted  using  this  formula. 

umx  =w(uW,  u$x  ,u£l,a)  umy  =  w(u$y,u$y,u$y,a)  (37) 

where  the  weighting  function,  uj  as  defined  by  Wang  and  Chang 


uj  = 


(■ 0m20m2.)aU'mxi  +  (0  m3#  ml)"  Umx  j  +  {0  ml®  U'mx , 

(0ml0m2)a  +  (9m2dm3)a  +  (0m30ml)a 


(38) 


where  Xi  =  x,y  and 


9 


mr  — 


(39) 


Equations  (35)  and  (36)  are  used  to  progress  the  odd  derivatives  to  the  next  time  step.  Since  the  odd 
derivatives  are  dependent  on  the  next  lower  derivative  it  must  be  calculated  before  hand. 


III.  Preliminary  Results 

To  date  we  have  simulated  the  convection  and  linear  acoustic  equations.  In  both  cases  a  4thorder  Taylor 
series  was  used  which  should  lead  to  a  convergence  rate  of  4.  To  determine  the  convergence  rates  we  used 
the  I2  defined  as 
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where 


£  —  V* numerical  ^ analytical •  (41) 

Since  the  Taylor  series  coefficients  are  constant  Eq.  40  is  simplified  to 

1 2  = 

where  Aj  is  the  area  of  cell  i. 

Each  convergence  test  was  done  on  two  different  types  of  meshes.  The  first  is  a  uniform  mesh  that  is 
created  by  splitting  a  quad  mesh  into  4  triangles.  The  triangles  are  created  by  splitting  the  quad  through  the 
diagonals.  The  second  mesh  is  an  unstructured  mesh  generated  using  an  advancing  front  algorithm.  Both 
meshes  were  created  using  Cubit.  For  the  uniform  mesh  the  QTri  mesh  was  used.  For  the  nonuniform  mesh 
the  TriAdvance  was  used.  The  TriAdvance  algorithm  first  tires  the  Delaunay  algorithm,  then  the  TriAdvance 
then  finally  the  QTri  method. 


III. A.  Advection  Equation 

The  governing  equation  is 

du  du  du 

TTT  +  - 1"  ay~^~  =  0. 

dt  ox  oy 

Under  periodic  boundary  conditions  and  a  domain  of  —7 r  <  x  <  7r,  —  7r  <  y  <  tt  u  is  equal  to 

u  =  sin(axx  +  avy  +  atf), 

where  ax  =  ay  =  1  and  at  =  —  yj a x  +  a^.  In  this  test  case  we  let  the  simulation  run  for  a  time  of  15  at  an 
approximate  CFL  number  of  1.0. 

When  calculating  the  convergence  rates  the  characteristic  length  was  the  square  root  of  the  average  area 
triangle  and  the  square  root  of  the  largest  area  triangle  for  the  uniform  and  nonuniform  mesh  respectively. 


cells 

h 

h 

0(h) 

400 

9.870E-02 

2.158E-01 

- 

1600 

2.467E-02 

6.726E-03 

5.00 

3600 

1.097E-02 

8.762E-04 

5.03 

6400 

6.169E-03 

2.065E-04 

5.02 

10000 

3.948E-03 

6.735E-05 

5.02 

(a)  Structured  grid 

Table  2:  Convergecne  rates 


cells 

h 

h 

0(h) 

232 

1.690E-01 

2.929E-01 

- 

926 

4.255E-02 

9.727E-03 

4.94 

2068 

1.909E-02 

1.315E-03 

4.99 

3694 

1.068E-02 

3.131E-04 

4.94 

5774 

6.835E-03 

1.022E-04 

5.01 

(b)  Unstructured  grid 

for  the  advection  equation. 


As  seen  in  Table  2  the  convergence  rates  are  5thorder  instead  of  4thorder  as  expected, 
investigating  this  in  more  detail  before  the  final  paper  is  submitted. 


We  plan  on 


III.B.  Linear  Acoustic 

The  governing  equation  is 


dp 

dt 


+  Po 


dvx\ 

dx  ) 


=  0 


dvx  al  dp 
dt  po  dx 
dvy  ah_dp_ 

dt  po  dy 
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where  p,  vx,  vy,  ao,  and  po  are  respectively  the  density,  x  and  y  component  of  the  velocity,  free  stream 
speed  of  sound  and  free  stream  density.  Under  periodic  boundary  conditions,  and  a  domain  —  2n  <  x  <  27 r, 
— 27t  <  y  <  27t  an  analytical  solution  is 


p  =  p'  sin(axa;  +  ayy  +  att) 

Vx  =  —a%  —  —wa.(axx  +  avy  +  att ) 

Po  at 

Vy  =  -al~—  sin(axx  +  avy  +  att ), 

Po  at 

where  p'  is  the  amplitude  of  the  perturbation.  In  the  test  case  we  set  ax  =  ay  =  1,  p'  =  0.01,  cto  =  rhoo  =  1.0. 
We  let  the  case  run  for  a  time  of  18  at  an  approximate  CFL  number  of  1.0. 

The  method  used  to  find  the  characteristic  length  is  that  same  as  the  one  used  for  the  advection  equation. 
As  seen  in  Table  3  the  convergence  rates  are  also  higher  than  expected.  It  should  also  be  noted  that  for 


cells 

h 

h 

0(Z2) 

cells 

h 

h 

0(h) 

1600 

9.870E-02 

4.189E-03 

- 

926 

1.705E-01 

4.330E-03 

- 

3600 

4.386E-02 

6.585E-04 

4.56 

2068 

7.636E-02 

7.193E-04 

4.47 

6400 

2.467E-02 

1.621E-04 

4.87 

3694 

4.275E-02 

1.756E-04 

4.86 

10000 

1.579E-02 

5.361E-05 

4.96 

5774 

2.735E-02 

5.849E-05 

4.92 

(a)  Structured  grid 

(b)  Unstructured  grid 

Table  3:  Convergecne  rates  for  the  acoustic  equation. 


similar  values  of  h  the  unstructured  mesh  tends  to  have  a  smaller  error  than  the  structured  mesh. 

IV.  Conclusion 

In  this  extended  abstract  we  presented  a  derivation  for  a  2D  4thorder  solver  on  an  unstructured  mesh 
using  the  CESE  scheme.  We  demonstrated  higher  order  convergence  using  the  advection  and  linear  acoustic 
equations. 

In  the  final  paper  we  will  also  include  solutions  of  the  Euler  Equation,  provide  a  computational  efficiency 
report  and  provide  an  explanation  as  to  why  we  were  able  to  achieve  a  convergence  rate  higher  than  expected. 
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LB  H 


M 


Fig.  1:  The  CESE  domain  using  a  triangular  mesh  with  points  of  interest  marked  by  squares,  circles  and 
crosses. 


(a)  Conservation  Element 


(b)  Solution  Element 


Fig.  2:  The  CE  and  SE  centered  on  the  centroid  mesh  point  (j,n). 
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V.  Figures 


( x0ly0ltn  1/2)r 


jr 

Fig.  3:  The  side  faces  associated  with  a  solution  point 
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4 


3 


Fig.  4:  The  bottom  face 


a 


f 


Fig.  5:  The  stencils  used  for  the  odd  derivatives 
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